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Abstract 

The  GPS  Master  Control  Station  (MCS)  currently  makes  significant  use  of  the  Allan  Variance. 
This  two-sample  variance  equation  has  proven  excellent  as  a  handy,  understandable  tool,  both  for 
time  domain  analysis  of  GI^  Cesium  frequency  standards,  and  for  fine  tuning  the  MCS's  state 
estimation  of  these  atomic  clocks. 

The  Allan  Variance  does  not  explicitly  converge  for  the  noise  types  of  a  < -3,  and  can  be  greatly 
affected  by  frequency  drift.  Because  GPS  Rubidium  frequency  standards  exhibit  non-trivial  aging 
and  aging  noise  characteristics,  the  basic  Allan  Variance  analysis  must  be  augmented  in  order  to  a) 
condensate for  a  dynamic frequency  drift,  and  b)  characteriie  two  additional  noise  types,  specifically 
a --3  and  a  = -4.  As  the  GPS  program  progresses,  we  will  utilize  a  larger  percentage  of  Rubidium 
frequency  standards  than  ever  before.  Hence,  GPS  Rubidium  clock  characterization  will  require 
more  attention  than  ever  before. 

The  three-sample  variance,  commonly  referred  to  as  a  renormalized  Hadamard  Variance,  is 
unaffected  by  Unear  frequency  drift,  converges  for  a  >  -S,  and  thus  has  utility  for  modeling  noise  in 
GPS  Rubidium  frequency  standards.  This  papa-  demonstrates  the  potential  of  Hadamard  Variance 
analysis  in  GPS  operations,  and  presents  an  equation  that  relates  the  Hadamard  Variance  to  the 
MCS*s  Kalman  Filter  process  noises  (qs). 


INTRODUCTION 

The  two-sample  variance,  or  what  we  commonly  refer  to  as  the  Allan  Variance,  has  been  an  excellent 
device  for  time  domain  characterization  of  GPS  Cesium  frequency  standards  over  the  past  few  years.  Over 
the  past  year,  the  GPS  Master  Control  Station  (MCS)  has  also  applied  the  Allan  Variance  towards  fine 
tuning  the  MCS’s  state  estimation  of  these  Cesium  clocks  [3]. 

In  terms  of  Power-Law  Spectral  Density  exponents,  the  Allan  Variance  does  not  explicitly  converge  for 
noise  types  of  a  <  -3,  and  may  be  greatly  affected  by  frequency  drift  [5].  Because  GPS  Rubidium 
frequency  standards  exhibit  significant  aging  and  aging  noise  characteristics,  the  Allan  Variance  analysis 
must  be  augmented  to  dynamically  compensate  for  this  frequency  drift,  and  to  characterize  two  additional 
noise  types,  specifically  a  =  -3  and  or  -  -4.  As  the  GPS  program  progresses,  we  will  utilize  a  larger 
percentage  of  Rubidium  frequency  standards  than  ever  before.  In  particular,  the  Block  HR  satellite 
platform  will  house  three  atomic  frequency  standards,  and  two  of  these  three  will  be  Rubidium.  Clearly, 
the  characterization  of  GPS  Rubidium  clocks  will  soon  require  more  attention  than  ever  before. 
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In  contrast,  the  three-sample  variance,  conunonly  referred  to  as  a  renormalized  Hadamard  Variance,  is 
unaffected  by  linear  frequency  drift,  converges  for  or  >  -5  [8],  and  hence  has  a  potential  utility  for  modeling 
the  various  noise  types  resident  in  GPS  Rubidium  frequency  standards.  This  paper  demonstrates  this 
potential  for  Hadamard  Variance  analysis  in  GPS  analysis  operations,  and  presents  the  relationship 
between  the  Hadamard  Variance  and  the  MCS’s  Kalman  Filter  process  noises  (^s). 


THE  HADAMARD  VARIANCE  EQUATION 


A  mainstay  of  atomic  clock  characterization,  the  two-sample  (Allan)  Variance  essentially  examines  the 
second  difference  of  phase,  equivalent  to  the  first  difference  of  the  time-averted  frequencies  over  two 
successive  adjacent  time  intervals  (r)  [5]: 


<r%(T)  = 


1 


2(M-i)tr 


yiCFivi  ~Pi)  >  Fj—the  time  -  averaged  frequency  over 


(1) 


Similar  in  principle  to  the  structure  to  the  Allan  Variance,  the  three-saavp\s  variance  examines  the  third 
difference  in  phase,  equivalent  to  the  second  difference  of  the  time-averaged  frequencies  over  three 
successive  adjacent  time  intervals  (r).  The  timing  community  has  commonly  referred  to  this  three-sample 
variance  as  the  Hadamard  Variance.  Though  the  term  Hadamard  Variance  has  been  used  more  generally 
in  various  applications  of  multi-sample  time  domain  analysis,  for  the  purposes  of  this  paper,  we  shall 
define  the  Hadamard  Variance  as  follows: 


H' 


r%(r)  = 


1 


M-l 


,  ><,.=016  time  -  averaged  frequency  over  r,. 


(2) 


The  Hadamard  Deviation  (the  square  root  of  the  Hadamard  Variance)  identifies  two  noise  types  that  the 
Allan  Deviation  does  not  explicitly  identify  [8],  For  this  paper,  we  shall  name  the  following  noise  types: 
for  a  =  -3,  “Flicker  Walk  FM”;  for  a  =  -4,  “Random  Run  FM”  [5],  Figure  1  visually  describes  the  noise 
types  identified  by  the  Hadamard  Deviation  [6], 


In  terms  of  phase,  equation  (2)  converts  to  (8]  ; 
1 


V-3 


y(r)  = 


6r^(W-3)tr 


y](Xj+3-3xj^2 +  3x,.^, -X,.)  ,  x,  =the  phase  measurement  at  f,. 


(3) 


or,  equivalently; 


^  cr^  (r)  =  ~  E[(X(^3  -  X(+2 )  -  (x..^^  -  x^+t )  -  (x,^2  ~x^\)  +  )f  (4) 

or 

where  E[.)  is  the  expectation  operator.  Each  phase  measurement  x;  =  x{ti)  in  equation  (4)  is  separated  fix>m 
each  neighboring  successive  phase  measurement  by  a  time  interval  value  of  r.  Meaning, 
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x{ti^,)  =  x{ti+  T), 


(5) 


xiti+i)  =x(ti+i  +  t)  =x(ff  +  2r),  and 


x(r,+j)  =  x(ti^i  +  r)  =  +  2it)  =  *(^,  +  3 1) 


(6) 

(7) 


MCS  KALMAN  FILTER  TIME  UPDATE  PREDICTIONS 

The  propagation  (time  update),  of  Rubidium  clock  states  in  the  MCS  Kalman  Filter,  is  modeled  using  the 
following  polynomial  expansion  [7]: 
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1 

T 

(l/2)r^ 

'x(t)' 

Ax 

yit+T) 

= 

0 

1 

T 

y(f) 

+ 

Ay 

zit  +  r) 

|0 

0 

1 

z{t) 

1 

Az 

(8) 


where  t  is  the  prediction  span,  and  x(r),  y(f),  and  z(f)  are  the  phase,  frequency,  and  frequency  drift  values, 
respectively,  of  the  clock  in  question.  Note  that  y(f)  is  the  time  derivative  of  x(t),  and  z(t)  is  the  time 
derivative  of  y(0-  A(x),  A(y),  and  A(z)  are  assumed  to  be  random  error  increments,  independent  of  x(r),  y(r), 
and  2(0,  having  a  prediction  covariance  P  represented  by  a  function  of  the  Kalman  Filter  process  noises 
(9s)[l,7]; 
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An  expansion  of  equation  (8)  produces  the  following  equations; 
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Using  equations  (10),  (1 1 ),  and  (12),  and  examining  the  differences  between  each  successive  x,: 
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(13) 


(*,+/  -  Xj)  =  +  {l/2)T^Zj  +  Arj+I 

(Xi+2  -  XfH./)  =  ty,+,  +  (l/2)/zi+;  +  Ax, >5  (14) 

(x,+j  -  Xi+j)  =  +  (1/2)t^z,+2  +  AXi+j  (15) 

AN  EXPECTATION  OPERATOR  EXPANSION  OF  fi(^y{z) 

Inserting  equations  (13),  (14),  and  (15)  into  the  following  expression: 

[ThirdDiJ]  =  [(x,+5-x,+2)  -  (x,+2-x*+0  -  (Xi^2-Xi^,)  +  (X(+;  -xj]  (16) 

obtains'. 

[ThirdDiJ]  =  Axj+j  -  2AXf+2  +  Ax,+/  +  r{0;  -3^^+;)  -  (>,+;  ->',>2)}  +  (l/2)z^{(Zi  - z,+,)  -  (Z(+;  -z,+2)}  (17) 

Examining  the  differences  between  each  successive >»,  and  Z(,  from  equations  (10),  (11),  and  (12): 

(yi+Z -Ji)^  s?r  + A>/,+;  (18) 

0.>2 ->',+;)  =  tri+/  +  Ayj+?  (19) 

(z,+,  -  z/)  =  Az/+/  (20) 

(^1+2  '  Zi+y)  =  Az,+2  (21) 

Equation  (17)  translates  into: 

[ThirdDiJ]  -  Ax, +2  -  2Ax,+2  +  Ax,+2  +  r{(i2i+/  +  A3',+2)  -  (tZ/  +  Ay,+;)}  +  (l/2)r^{(Az,+2)  -  (Az,>;)}  (22) 

With  some  more  algebraic  manipulation: 

[ThirdDiJ]  =  Ax(+j  -  2AXj+2  +  Ax,+;  +  r{(A3;,+2  -  Ay,+;)  +  (jZi+/  -  zr*)}  +  (l/2)r^{(Az,+2)  -  (Azi+/)}  (23) 

[ThirdDiJ]  =  Ax,+2  -  2Ax,>2  +  AXi+/  +  r{(A3>,+2  •  Ayj+2)  +  (rAzi+;)}  +  (l/2)t^((AZ(+2)  -  (Azi+;)}  (24) 

[ThirdDiJ]  =  Ax, +2  -  2Ax,+2  +  Ax,+;  +  r{(Ay,+2  -  ^3-,+;)}  +  (l/2)z^{(Az,+2)  +  (Az,+;)}  (25) 

[ThirdDiJ]  =  (Ax, *2}  +  {-2Ax,.2  +  t^y^^2  +  (l/2)z^(Az,+2))  +  (Ax,.,  -  rAy,.,  +  (l/2)t^(Azi.,)}  (26) 

Since  we’ve  now  broken  down  this  expansion  into  three  independent  polynomial  terms: 

E[rhirdDjyi^  =  E[( Ax,+2}  +  {-2Ax,+2  +  zAy,+2  +  (l/2)7^(A2,+2)}  +  {Ax,+,  -  zAy,.,  +  (1/2)t^(Az,./)}]*  (27) 

the  independence  of  each  term  {}  allows  us  to  separate  equation  (27)  into  three  individual  expectation 
operators  [4]: 
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WrhirdDifi^  =  E[AX(^5]"  +  E[-2AXi,,  +  +  (1/2) ^Az, >2)]"  +  E[Ax,+,  -  rAy,>;  +  (l/2)i2(Azi.,)]'  (28) 

Expressing  each  term  ofU[ThirdDiJ^  as  a  function  of  the  Kalman  Filter  prediction  covariance  matrix  [IJ; 


EfAr,,3f  =[1  0  0]  2^ 


=  ?,  r  +  (1  /  3)^2  +•  (1  /  20)^3  r* 


E[-2Ax,,,+rA>-,,2+(l/2)r^(Az,,2)f  =[-2  r  (l/2)r^]  P  r 

[(1/2)^' 

=  49,  r + (1  /  3)^2  +  (9  /  20)^3  r  * 


E[Ax„,-rAy,,,+(l/2)r'(Az,,,)f  =[l  -r  (1/2)t^]  P  -r 

L(l/2)r' 

=  g,r  +  (1  /  3)^2  +  (1  /  20)^37^ 


(29) 


(30) 


(31) 


By  adding  each  term: 

EVThirdDijf  =  69,(7)  +  92(1^)  +  (1  l/20)93(r^)  (32) 

Hence, 

;,cr%(r)  =  ^E[77»>*rfD//f  =^[69,(r)  +  92(^')  +  (ll/20)93(r^)]  (for  a=0,  -2.  -4)  [6]  (33) 

or  or 

„a^^r)  =  9i  f '  +  (1/6)927  +  (1 1/120)93:^  (for  a  =  0,  -2,  -4)  (34) 


RELATING  WHITE  PM  TO  THE  HADAMARD  VARIANCE 

Equation  (34)  does  not,  however,  account  for  white  PM  noise  (a  =  2)  [6],  sometimes  also  referred  to  as 
representation  error  [1],  The  Hadamard  Variance  can  be  expressed,  in  terms  of  phase  measurements,  as 
follows: 

H =  -!yE[x,.^j  - 3x,.^2  +  " Xif  (35) 

67 

When  white  PM  is  the  only  significant  noise  component,  the  individual  phase  values  are  uncorrelated  with 
time,  and  may  be  separated  [4]; 
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^  cr%  ( r )  -  ^  {E[X(^3  ]^  +  E[-3Jf,+2  f  +  E[3x,^,  ]^  +  E[-X(  ]^ }  (for  a  =  2) 


(36) 


When  assuming  that  white  PM  is  the  primary  noise  source,  the  representation  error,  which  we’ll  denote  as 
qo  ~  E[x,]^  is  independent  of  and  thus  is  a  phase  variance  that  is  constant  across  time.  Therefore, 


„flr%(r)  =  ~{(l4-9-f-9+l)E[x,f}  =  ^?o  (fora=2) 


(37) 


THE  HADAMARD-Q  EQUATION 

In  the  presence  of  both  a)  white  PM,  and  b)  the  three  noise  types  modeled  by  P,  and  assuming  independence 
between  the  white  PM  and  the  other  noise  types,  equations  (34)  and  (37)  can  be  combined  into  one  that 
models  four  noise  types,  namely  a  =  2,0,  -2,  and  -4; 

=  (10/3)9of'  +  qif'  +  (1/6)^2T+  (1 1/120)93 (38) 

Note  how  this  equation  compares  to  the  analogous  equation  relating  the  Allan  Variance  to  the  9s  [2,3]; 

=  2qof^  +  9if'  +  {1/3)92T+  (1/20)93^^  (39) 

For  white  FM,  the  Allan  and  Hadamard  Variances  are  mathematically  equivalent.  For  white  PM,  the  two 
Variances  are  roughly  the  same,  and,  for  random  walk  FM,  the  Variances  differ  by  a  factor  of  two. 
Though  an  analyst  may  use  either  the  Hadamard  Variance  or  the  Allan  Variance  for  deriving  MCS  qs,  each 
has  its  own  set  of  advantages  and  disadvantages. 

The  primary  advantage  of  the  Hadamard  Variance  is  the  automatic  removal  of  linear  frequency  drift  [8]. 
Whereas  the  equation  relating  the  Allan  Variance  to  MCS  9s  assumes  that  the  analyst  must  apply  a 
continuously  dynamic  correction  for  frequency  drift,  the  Hadamard-Q  equation  doesn’t  require  this 
assumption.  The  tradeoff,  however,  is  that  the  Hadamard  Variance  incurs  an  extra  computational  burden, 
simply  because  it  examines  the  third  (vice  the  second)  difference  of  phase.  For  analyzing  GPS  Cesium 
frequency  standards,  the  increased  computational  load  of  the  Hadamard  Variance  proves  fruitless,  only 
because  the  Allan  Variance  gets  the  job  done  more  efficiently  [3]. 

Many  timing  experts,  over  the  years,  have  extensively  used  techniques  for  applying  a  continuously  dynamic 
correction  for  frequency  drift,  prior  to  using  the  Allan  Variance  for  deriving  93  values  with  high  confidence. 
This  paper  does  not  address  the  issue  of  confidence  in  the  93  value  produced  by  frie  Hadamard-Q  equation. 
This  paper  does,  however,  present  an  easily  understood  relationship  between  a  relatively  lesser  known 
equation  (the  Hadamard  Variance),  and  a  set  of  system  parameters  used  by  the  MCS  (the  Kalman  Filter 
qs).  On  initial  appearance,  given  the  computational  capability,  one  can  see  the  great  potential  utility  of  an 
algorithm  that  applies  a  relatively  simple  equation  onto  a  large  measurement  data  base,  in  order  to  derive 
Kalman  Filter  9s,  without  the  need  to  apply  preparatory  frequency  drift  corrections.  In  the  future,  the 
author  hopes  to  further  investigate  a)  the  real-world  utility  of  this  relationship,  b)  the  issue  of  estimate 
confidence,  and  c)  the  net  gain  from  the  increased  utility  balanced  against  the  increased  computational 
burden. 
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CONCLUSION 


The  implication  of  the  Hadamard  Variance  in  GPS  operations  is  as  follows:  Analysts  at  the  MCS  could 
simply  gather  a  large  data  base  of  clock  phase  measurements,  apply  all  known  step  corrections,  perform  a 
number  of  iterations  of  the  Hadamard  Variance  equation,  and  plot  the  results  to  visually  describe  the  noise 
characteristics  of  GPS  Rubidium  clocks  (including  a  =  -4).  Consequently,  the  Hadamard  Variance  could 
offer  GPS  operators  an  alternate  tool  for  characterizing  GPS  atomic  frequency  standard  noise. 

Perhaps  more  significantly,  the  implication  of  equation  (38)  is  that  GPS  analysts  now  have  an  easily 
understood  technique  to  relate  raw  clock  phase  measurements  towards  deriving  important  Kalman  Filter 
clock  estimation  parameters  (^s),  that  are  unique  to  the  performance  of  each  individual  clock,  and  that  will 
include  qs  automatically,  without  any  need  for  the  preliminary  removal  of  frequency  drift. 

The  MCS  hopes  to  make  continued  use  of  the  Allan  Variance,  for  both  the  characterization  of  Cesium 
clocks,  and  the  derivation  of  their  associated  process  noise  values.  This  application  of  the  Hadamard 
Variance  widens  the  array  of  available  tools  for  the  characterization  of  Rubidium  clocks,  and  the  derivation 
of  their  associated  process  noise  values.  The  GPS  Block  HR  satellite  program  will  use  a  large  percentage 
of  Rubidium  clocks.  Certainly,  the  ever-important  issue  of  refining  Rubidium  clock  estimation  may  see  its 
most  important  days  in  the  years  ahead. 
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Noise  Types  Identified  by  the  Hadamard  Deviation 


Noise  Type _ Power>Law  Spectral  Density  Exponent  HCTy(T)  Log-Log  Slope 


White  PM 

a  =  2 

-1 

Flicker  PM 

a=  1 

-1 

White  FM 

a  =  0 

-1/2 

Flicker  FM 

a  =  -l 

0 

Random  Walk  FM 

a  =  -2 

1/2 

Flicker  Walk  FM 

a  =  -3 

1 

Random  Run  FM 

a  =  -4 

3/2 

Figure  1 
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Questions  and  Answers 

SERGEY  V.  ERMOLIN  (HEWLETT-PACKARD):  The  Hadamard  Variance  does  remove 
linear  drift,  that’s  true;  and  that  saves  you  some  time  on  preprocessing.  But,  it  doesn’t  remove 
any  drift  beyond  linear.  As  you  showed  on  one  of  your  first  slides,  that  rubidium  standards  on 
board  one  of  the  space  vehicles  show  not  only  linear  drift,  but  possibly  quadratic  drift  to  some 
other  power. 

So  still,  if  you  wanted  to  go  ahead  with  longer  operational  time,  you  would  still  have  to  do 
some  preprocessing,  even  if  you  used  the  Hadamard  Variance. 

CAPTAIN  STEVEN  HUTSELL  (USAF):  The  analogy  is  the  Allan  Variance  does  not 
care  about  a  constant  frequency  offset.  For  instance,  our  atomic  clocks  can  have  a  1  x  10“^^ 
frequency  offset;  but  if  it’s  stable  enough,  the  Allan  Variance  will  be  low,  regardless  of  that 
frequency  offset.  The  analogy  is  in  the  Hadamard  Variance,  if  the  satellite  clock,  or  whatever, 
has  an  already  existing  common  offset  of,  say,  3x  10~^®  seconds  per  second  squared  of  frequency 
drift,  that  will  not  adversely  affect  the  Hadamard  Variance  calculations. 

However,  in  the  Allan  Variance,  if  there  is  a  frequency  drift,  it  will  affect  it.  But  with  the 
Hadamard,  it  won’t.  In  the  same  sense,  a  random  walk  in  random  walk  FM  will  affect  the 
Allan  Variance  by  causing  a  positive  one  slope.  But  the  frequency  offset  itself  will  not  affect  it. 

DAVID  ALLAN  (ALLAN’S  TIME):  Dr.  Barnes  did  some  work  three  decades  ago  on  the 
confidence  question,  actually  related  to  it,  in  the  third  difference  estimate.  The  bottom  line 
was  that  the  confidence  is  worse  by  a  significant  amount,  especially  when  you  have  finite  data 
lengths  that  come  into  impact  you  quite  adversely. 

The  other  point  is  that  it  turns  out  a  logarithmic  drift  estimator,  both  for  quartz  and  rubidium 
is  much  better  than  linear.  I  think  the  graph  that  you  showed  outside  of  the  turn-on  transient 
probably  would  fit  a  logarithmic  curve  quite  well. 

So,  one  might  be  better  doing  logarithmic  modeling  if  we  deal  with  a  lot  of  rubidiums  in  the 
future.  Yes,  one  would  expect  that  the  logarithmic  function  could  be  fit  to  this  quite  well, 
outside  of  the  first  point. 

The  last  one  I  would  like  to  suggest  that  you  think  about  —  and  I  want  to  say  that  I  think 
you’ve  done  a  beautiful  piece  of  work,  but  you  can  actually  remove  the  effect  of  drift,  kind  of 
in  real  time,  from  the  second  difference  operator,  because  you  know  the  exact  equation  for 
the  effect  of  drift  on  it  until  you  can  subtract  that  from  and  get  an  estimator  variance  without 
the  drift  effectively  in  real  time.  So,  it  doesn’t  need  to  impact  the  value  of  the  variance  if  you 
don’t  want  it  to;  and  it  gives  you  a  tighter  confidence  of  the  estimate. 

CAPTAIN  STEVEN  HUTSELL  (USAF):  Yes,  and  the  intent  of  this  is  not  really  to  present 
the  best  way  to  estimate  frequency  drift  dynamically.  Really,  what  I  v/anted  to  do  was  examine 
the  way  the  MCS  is  currently  set  up.  And  right  now,  it  does  not  have  the  capability  to  do 
what  you  just  described.  The  MCS  is  only  set  up  as  a  dynamic  Kalman  filter  three-state  vector 
that  needs  process  noise  values. 

I  completely  agree  that  are  far  more  sophisticated  techniques  to  look  at  frequency  drift  than 
what’s  set  up  in  the  Kalman  filter.  Sorry,  your  first  comment? 
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DAVID  ALLAN  (ALLAN’S  TIME):  [Inaudible]. 

CAPTAIN  STEVEN  HUTSELL  (USAF):  Yes,  we  see  that  too.  We  see  it  start  to  converge 


DAVID  ALLAN  (ALLAN’S  TIME);  [Inaudible]. 

CAPTAIN  STEVEN  HUTSELL  (USAF):  We  agree.  Over  time,  the  frequency  drift  goes 
from  a  negative  value,  around  3  or  4  x  10“^®,  and  gradually  starts  logarithmically  to  approach 
zero. 

At  the  beginning,  however,  sometimes  we  see  it  start  hugely  negative,  like  -lx  Sometimes 

-f-1  X  We’re  talking  about  over  the  first  48  hours  that  we  turn  it  on.  We  would  need  to 

address  how  we  try  to  model  that.  It’s  also  probably  appropriate  to  ask  what’s  causing  that, 
what  physical  phenomenon  is  causing  it  to  be  positive  at  the  beginning  for  some  clocks,  and 
negative  for  the  others.  But,  I  do  agree. 
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